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Abstract 

We study the distribution of the largest eigenvalue in formal Hermitian one-matrix 
models at multicriticality, where the spectral density acquires an extra number of 
k — 1 zeros at the edge. The distributions are directly expressed through the norms of 
orthogonal polynomials on a semi-infinite interval, as an alternative to using Fredholm 
determinants. They satisfy non-linear recurrence relations which we show form a Lax 
pair, making contact to the string literature in the early 1990's. The technique of 
pseudo-differential operators allows us to give compact expressions for the logarithm 
of the gap probability in terms of the Painleve XXXIV hierarchy. These are the higher 
order analogues of the Tracy- Widom distribution which has k = 1. Using known 
Backlund transformations we show how to simplify earlier equivalent results that are 
derived from Fredholm determinant theory, valid for even k in terms of the Painleve II 
hierarchy. 



1 Introduction 



Random matrices appear in many branches of mathematics, physics and other sciences. 
Perhaps the best example to illustrate such a breadth of applications is a single formula 
known as the Tracy- Widom (TW) distribution, and we refer to [1] for an incomplete list. 
It was first realised as the scaled distribution of the largest eigenvalue of a large complex 
Hermitian random matrix with Gaussian entries [2] that can be expressed as the derivative 
of the Fredholm determinant of the Airy-kernel [3]. Apart from its appearance in many 
other areas this distribution is remarkably robust in the sense that it holds for a much larger 
class of non-Gaussian distributions, of either invariant [I] or Wigner type [5]; this property 
is known as universality. 

Another striking feature is the relation between the TW distribution and integrable 
hierarchies of non-linear differential equations, notably the Painleve hierarchies. A relation 
between matrix models and integrable hierarchies was first observed in the string literature 
of the early 1990's where matrix models were applied to two-dimensional quantum gravity, 
and we refer to [B] for an extensive review. There the possibility of having multicritical points 
was introduced in order to couple matter to gravity. This happens when k — 1 additional 
zeros collide with the spectral density - which is the semi-circle in the Gaussian case (k = 1). 
Such scenarios can be realised by fine-tuning non-Gaussian potentials, see e.g. [7J |8] for a 
classification and references. In order to have the largest possible set of critical exponents 
available this also included so-called formal matrix models, where the confining potential is 
not bounded from below. This allows k to assume any integer value. 

Whereas the limiting kernel of orthogonal polynomials with multicritical potentials has 
been considered before in mathematics, e.g. in [HJ [TO], and in the physics literature, see [TT] 
for a comprehensive list of references, it was not until recently that the question about the 
existence of higher order analogues of TW distributions at such multicritical points was first 
answered |12J . This also includes a first numerical investigation in [13]. Building on the 
Riemann-Hilbert analysis of the limiting kernel for k = 3 [9], in [12] Fredholm determinant 
theory was used to describe the higher order analogues of TW for any odd k 0. They were 
found to be characterised by a system of equations including the odd members of the Painleve 
II hierarchy. 

Several alternative derivations of the TW distributions exist, in particular a recent heuris- 
tic but very transparent one [H], that avoids the Fredholm determinant formalism. Starting 
from a Gaussian potential they directly compute the gap probability that the interval {y,oo) 
is empty of eigenvalues. It is given in terms of the norms of orthogonal polynomials on the 
complement (— oo,y] which satisfy a set of non-linear differential equations. After carefully 
performing the asymptotic analysis of the recursion coefficients in the three-step recurrence 
relation TW at k — 1 follows. 

Our goal is to generalise the approach |14j to the higher order multicritical cases. How- 
ever, it rapidly becomes very cumbersome, already when extending their derivation of TW 
to a generic quartic potential (which only shows its known universality). The main idea in 

^^Note that our convention of counting zeros differs from |12j . 
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this paper is to cast the recurrence relations into a different set of flow or string equations, 
that allow to make close contact to integrable hierarchies via the Lax pair. Once these com- 
mutation relations become multiplicative we can apply the Lax method to solve them using 
pseudo-differential operators, following [61 [15] . Ultimately this leads to a hierarchy of differ- 
ential equations given by the Painleve XXXIV hierarchy for any integer k, in terms of which 
the logarithm of the gap probability at multicriticality can be expressed. Here we exploit 
that the limiting case y — > oo is well understood in terms of these integrable hierarchies, see 
e.g. [6j. Furthermore, we manage to prove the equivalence to the subset of odd- A; hierarchies 
that were derived rigorously in [12J. The link is provided by a Backlund transformation that 
was explicitly constructed in [IB] . Applied to our result it leads to all members of Painleve 
II including [12J. A shift property of the so-called Lenard differential operators (also known 
as Gelfand-Dikii polynomials) allows us to simply their result. 

The remainder of this paper is organised as follows. In section [2] we review the orthogonal 
polynomial formalism adopted to [H] and derive our flow and string equation for finite matrix 
size. The main result is derived in section [3] where the double scaling limit is taken on the 
string equation, leading to Painleve hierarchies for the gap probability. Explicit examples 
including TW are given and the equivalence to (12] is shown, before concluding in section HI 
Various technical details are delegated to the appendices. 



2 The Orthogonal Polynomial Formalism 

Consider the following truncated partition function 

N N 



Z N (y;a,{ gi }) = -1 f f[ dX, e~ NaV ^ f[(X k - X 

Iy -J-oo i=1 k ■ 



2 

k>j 



(2.1) 



V W = Ey^- ( 2 - 2 ) 

1=1 

Here V(X) is a formal power series and a > is a real parameter. Both are independent 
of y. When sending the upper integration range y — > oo it corresponds to the standard 
Hermitian one-matrix model Zjv(oo; a, {<?;}), given in terms of eigenvalues A, of a random 
matrix M. In general we want to compute the probability P7v(A max < y; a, {gi}) that the 
largest eigenvalue of the random matrix M, is less than y. It has the following expression, 

to f \ / j u Z N (y;a,{gi}) 

IVA m ax < V, a, {gi\) = —— p-r-, (2.3) 

Z N (oo;a,{gi}) 

which is why our definition eq. (12. ip is convenient. One possibility to determine this prob- 
ability is to introduce orthogonal polynomials for the partition function Zjv(oo; a, {gi}) and 
to study their asymptotic behaviour in an appropriate scaling limit. However, we will follow 
[H] here by introducing polynomials for the truncated partition Zjsf(y] a, {gi}). Formally 
this could also be formulated in the one-matrix model by adding a hard wall at y to the 
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potential. Introduce a set of polynomials {7r n (A) :n6N} such that 7r ra (A) is of order n and 
they are orthonormal with respect to the inner product defined by, 

(7r n | n m ) = [ V d\ e- NaV ^7r n (X)7r m (X) = 5 nm . (2.4) 



Note that the coefficients in the expression for each ir n depend on y, a and the potential. 
Furthermore, we define h n = h n (y) > by the leading coefficients 

vr n (A) = -L X n + 0(A n_1 ) . (2.5) 

They are positive being the squared norms of the polynomials in monic normalisation. It 
can be shown using the orthonormality that p2] 

N-1 N-1 

Z N (y; a, {g<}) = ]J h = h o U ■ ( 2 -6) 

t=0 i=l 

Here we have defined the ratios 

V_ f or n > 1 , (2.7) 



h 



n-1 



that will determine the distribution of the largest eigenvalue. It is useful to rewrite 
Z N (y;a,{gi}) as, 

N-1 , . v 

teg[Z N (y ] a,{gi})]=N]og[h ] + N^(l-^Uog[r i ]. (2.8) 

i=l ^ ' 

In the case of a Gaussian potential gi = 25 2) i, and y — > oo everything is explicitly known for 
finite- N because the n n are proportional to Hermite polynomials, with h n = n\y/7r/2 n and 
the ratio assuming the simple form r n = n/2 (we have set Na = 1 for simplicity). 

In order to determine the r, in a general setting we first introduce a set of multiplication 
and differentiation operators on the ring of polynomials n n which form a complete set of 
functions. We begin with the multiplication operator B: 

B nm n m (\) = \n n (\) = y/r n+1 n n+1 (X) + s n 7r„(A) + y/r^iTn^X) , (2.9) 

where we have used summation conventions for the matrix multiplication on the left hand 
side. In terms of the inner product this reads 

The fact that B is tridiagonal follows from the three step recurrence relation that arbitrary 
orthogonal polynomials with weights on the real line satisfy; this includes the case of the 
truncated integral in ( 12. lft . Besides r n we encounter here a second recurrence coefficient, 
s n = (*7r n | A7r n ), that must be determined in principle. Obviously B is a symmetric matrix 
as B nm = (ir n \ \iT m ) = (7r m | A7r„), or in operator language B = B T . 
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Differentiation with respect to the argument of the polynomials is denoted by A: 

A nm TT m (X) = <9 A 7T n (A) , (2.11) 

or in terms of the inner product A nm = {^ m \dxTT n ). Because the coefficients h n are A- 
independent the polynomial on the right hand side is at most of degree n — 1, and thus 

A n , n+k = for k > -1 , (2.12) 

or in other words A is strictly lower triangular. In particular for the uppermost non-zero 
diagonal we have 

A n>n - X = I" d\ e - NaV ^n n ^d x (-L X n + 0(A"~ 1 )) = n/^T n . (2.13) 

J — oo V V "n J 

For the Gaussian potential other explicit expressions can be derived. It is easy to see that 
A and B form canonical commutation relations 

[B,A] = 1, (2.14) 

due to (BA — AB) nm Tr m = d\(\ir n ) — \d\n n = n m 5 ntm , upon using the definitions. 

Finally we also introduce differentiation with respect to the truncation y denoted by C: 

C nm -K m {\) = d y ic n (\) , (2.15) 

or equivalently C nm = (ir m \ d y n n ) . Now the leading order coefficient containing h n in eq. 
(12. 5 j) is y-dependent, and so the polynomial on the right hand side is at most of degree n, 

C n:n+k = for k>0. (2.16) 

The matrix C is lower triangular including the diagonal, which is given by 

C nn = d\e- NaV ^TT n d y (^= A n + 0(A n - 1 )^) =~d v ]og[h n ] ■ (2.17) 

For further explicit expressions for the matrix elements see eq. (I2.2ip below. Matrix C 
satisfies the following commutation relations with the multiplication operator B 

[B,C] = -d y B, (2.18) 

where we have been careful to remember the y-dependence of the recursion coefficients, 
[BC - CB) nm iT m = B n idy-ni - d y (X'K n ) = B rd d y Tii - dy(B n i-ni). We are now prepared to 
introduce two composed operators that will play a crucial role in the following, by establishing 
a link to integrable systems. 
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2.1 The Flow Equation 



We begin with the matrix P, 



aN ( aN \ 

Pnm = A nm + C nm - —V\B) nm = (7T m | f d X + d y - — V' (X) J 7T n ). (2.19) 

It is a simple matter of integration by parts to see that P is anti-symmetric, P + P T = 0: 
P nm = f V d\e- NaV ^ m (d x + d y -^V'(X)]n 



2 

e- JVay ^7r w (y)7r„(y) - f rfA e"^^ (-AW(A)7T w 7T n + n n d x n m ) 

J — oo 

dA e _,Va7(A) f d y (7r„7r m ) - 7T n ^7r m —V'(X)ir m TT r 



dX e- NaV ^7r n (d x + d y - ^V'(X) ) vr m = -P mn . (2.20) 



Here we have used the fact that 



= d y 5 nim = d y (7r n \7T m )= e- NaV ^7c n (y)7f m (y)+ dX e~ NaV ^d y (7r n 7r m ) 
^0 = e- NaV ^n n (y)7T m (y) + C mn + C nm (2.21) 



which we have also expressed in terms of the matrix elements of C . Secondly, using (I2.12p 

Na 
2 

symmetry of P we find that 



and (12.1 6p we have that P n .n+k = —^V'{B) n>n+ k for all k > 0, and hence using the anti- 



P = ~(y'(B) + -V'(B).) f (2.22) 

where + and — denote the upper and lower triangular parts of a matrix respectively. Ob- 
viously its diagonal part vanishes, and its antisymmetry (given the symmetry of B) is now 
manifest. Because we have expressed P only in terms of powers of B, which is tridiagonal, 
we see P has only a finite number of non-zero off diagonals. For the Gaussian potential with 
gi = 262,1 (that is V(X) = A 2 ) for example it reads 

Pnm I Gauss 

—Na(^r n+ i <5 n +i,m — aA™<^-i,»j) • (2.23) 

Finally, given the commutation relations (12. 14ft and (12.18!) we arr ive at 

[P, y -B] = d y (y-B) = l- d y B. (2.24) 

This is our first main result of this section, and we will refer to it as the flow equation. 
Note that this equation implies that the matrices y — B and P form a Lax pair, where we 
have added a trivial zero, = [P, y], to bring it to that form. The flow equation (12. 24ft 
provides explicit relations among the recurrence coefficients r n and s n , including derivatives 
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with respect to y. To give an example we again consider the Gaussian case. Inserting (I2.10p 
and P in the form of (12.231) into (I2.24p we obtain five equations. Those for m = n ± 2 are 
identically satisfied, the equations for m = n ± 1 are the same, and thus we are left with the 
following two equations for indices m = n — 1 and m = n: 

r n+l ~ r n 

The first of this set of equations already appeared in [14] . however the second is new. The 
equivalence of our flow equations to [14J is shown in Appendix A. The advantage of our 
approach is that we may obtain such recursion relations for any potential. Furthermore, we 
can see that our recursion relations enjoy an explicit link to integrable systems defined by 
the Lax pair y — B and P. Indeed we can see this connection in this particular example as 
it is trivial to eliminate s n , thereby obtaining, 

Vi + rw - 2r„ = ^J^d* log[r B ] . (2.27) 

By substituting r n = exp[0„] we see that this is equivalent to the difference of the Toda 
lattice equation at two neighbouring lattice sites. 



2aN 
2aN^ 



d y \og[r n ] , 



(2.25) 
(2.26) 



2.2 The String Equation 

The second operator we introduce is the matrix H. It will lead to a so-called string equation 
that is algebraic in r n and s n . It is defined as, 

Na , 1 
H nm = (A{B-y)) nm —(V'(B)(B - y)) nm + -5 nm 

= J" d\ e- NaV ^n m ((A - y) (d x - ^V'(X)) + n n . (2.28) 

It is easy to see that H is antisymmetric, H + H T = 0, by applying the following identity in 
terms of a total derivative: 

= f dXd x [(X-y)e' NaV ^\ n (X)rc m ^)} 

J — oo 

dX e -N«v(x) [nn7rm + (A _ y )(_ iV ay'(A)7r n 7r m + (9 A 7r„)vr m + n n (d X 7c m ))} 

Hnm H mn . (2.29) 

Because A is strictly lower triangular and B tridiagonal we have that its product does not 
contribute to the + part of H, hence H n ^ n+ k = —^y l (V'(B)(B — y)) n ^ n+ k for all k > 0. Due 
to the antisymmetry of H we thus arrive at 

H = -^({V'{B){B-y)) + -{V'{B){B-y))_) . (2.30) 



oo 

y 
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This shows that also H only has a finite number of non-zero off-diagonals, being a multiple 
of B in its + and — part. For the Gaussian potential it reads 

H nm | Gauss + y/ r n+l(s n +i + S n — y)5n+l,m) 

+Na( y /r n ~(s n + s n _i - y)5 n -x, m + v / VW^-2,m) • (2.31) 

Furthermore, from the definition (12. 28ft and the commutation relations ( 12.14ft and ( 12. 18ft we 
obtain 

[B-y,H] = B-y , (2.32) 

which is the second result of this section (where we have added a trivial zero, = [—y, H]). 
We shall refer to it as the string equation and it provides relations among the recurrence 
coefficients, for any potential V, which are algebraic and no longer contain derivatives. Let 
us give again the Gaussian example. Plugging eq. (12. 101) and (12.3 II) into the string equation 
we obtain five equations, two of which for m = n ± 2 are identically satisfied. The ones for 
m — n ± 1 are identical and we are left with the two equations for m = n + 1 and n = m 
which read, 

S n+1 ~ S n + r n+2 ~ T n ~ y{s n +l — Sn) 

r n+1 (s n+1 + s n -y) - r n (s n + s n ^ - y) 

Finally we note that the two operators P and H are related. Because multiplication of 
matrices and taking their upper (+) or lower (-) part in general do not commute we derive 
the following relation in Appendix [Bj 

H = ( B -y)P-{BP) d -{B + ~B_)^V>{B) d 

= \{B,P}-yP-^{B + -B_,V'(B) d } , (2.35) 

where the subscript d denotes the diagonal part of a matrix. In the second line we have 
made the anti-symmetry of H manifest, by decomposing it into its (vanishing) symmetric 
part and its anti-symmetric part. This relation will be important when taking the double 
scaling limit in the next section. 



1 



Na 1 
1 

2iV^ ( 



y) 



(2.33) 
(2.34) 



3 Double Scaling Limit and the Lax Method 

In this section we will take the large-iV limit and determine the asymptotic recurrence 
coefficients. This will give the distribution of the largest eigenvalue by differentiating (12. 3p . 
In the case of a Gaussian potential and y — > oo it is well known that the limiting spectral 
density p(X) is given by the Wigner semi-circle that vanishes as a square root at its end 
points. By taking the large- N limit while simultaneously zooming into the vicinity of the 
right end point and letting y approach the same end point - hence the name double scaling 
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limit (d.s.l) - the distribution of the largest eigenvalue can be found. In [2] this was done 
using in an approach using the truncated partition function (12. ip which we generalise here. 

For non-Gaussian potentials it is well known [7J |S] how to construct multicritical poten- 
tials Vk y k'i such that the spectral density behaves as p(A) ~ (2a + — \) k ~ 1 / 2 (\ — 2a_) fc '~ 1 / 2 . 
Here (2a_, 2a + ) is the support of p which we will assume to be of single interval or one-arc 
type in all the following, and the integers k, k' — 1, 2, . . . label the degree of multicriticality, 
following ref. [8|. So for the Gaussian potential we have (k, k') = (1, 1) with a + = a_ = 1 
0. For these critical points a non-trivial d.s.l can be taken if we make the following scaling 
ansatz [TH] , 

f = ic-a k x (3.1) 
N = a^ 2k+1 ^ 2 , 

where a is the scaling variable given in fractional powers of iV and n/N — > £ becomes a 
continuous variable □. In addition we make the following general ansatz for our recurrence 
coefficients r n and s n : 

r{£,,y) = r c (l + ap x {x, s) + a 3/2 p 2 (x, s) + . . .) , 

s{£,y) = y/r^isc + aa^x, s) + a 3/ V 2 (x, s) + . . .) . (3.2) 

It will also be useful to introduce the functions, 

Ui = pi + Oi , v i = pi - ai . (3.3) 

In [4] it was shown that for a one-arc support and analytic potentials the r n and s n indeed 
become single valued functions, and we assume in the following this remains true for finite y. 
The d.s.l in the case y = oo is fully reviewed in Appendix O However, to quickly recap here; 
by substituting the above scaling ansatz into the y = oo string equations found in Appendix 
ICj or alternatively using pseudo-differential operators to solve the same string equation, one 
finds the result for k' < k, 

C k [ Ul {x,oo)} = x/2, (3.4) 
C k '[vi{x, oo)] = 0. (3.5) 

The Ck are the Lenard differential operators introduced in the next subsection, see (I3.16p . 
In the remainder of this section we argue that the scaling ansatz given in (13. 2p together with 
(13. ip . augmented to, 

f = l-a k x 

y = yc + acxs (3.6) 
N = a-^ 2k+1 ^ 2 , 

where C\ is an arbitrary constant to be fixed later, and s is the microscopic variable that 
describes the distribution of the largest eigenvalue around the right endpoint of the support 

2 Note that in [T2] the Gaussian case is labelled by I = k — 1 = 0, and that only critical points with our k 
taking odd values are considered. 

3 Note that we always have £ c = 1. This will be justified later. 
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y c (= 2a+), gives a non-trivial scaling limit in the case of finite y. Furthermore for general 
k' < k we derive the equation satisfied by u\ and argue that the equation (13. 5p is still valid. 
Each of these results are also derived in Appendix [D] in a more rigorous manner. 

For a particular potential V^k' we could of course follow [H] and insert (13. 6p and (13. 2 j) 
into either the flow or string equation ( I2.24[) or fl 2 . 3 2 [) and then Taylor expand in order to 
derive the relevant differential equations for the scaling functions ttj and Vi. However for 
general k we require a general method that can uncover the structure of such a hierarchy of 
equations. The particular approach we adopt here is to follow the method of Appendix O 
and solve (I2.32j) using pseudo-differential operators. For a review of this technique we refer 
to [6J as well as to [15] for its relation to the Lax method. 

3.1 Scaled String Equation and Lenard Differential Operators 

In the large- N d.s.l our matrices B, P and H which all have only a finite number non-zero 
off diagonals will scale to differential operators of finite order. The commutation relations 
that we derived in the previous section will then allow us to directly obtain a differential 
equation for the scaling function. First consider the scaling limit of the matrix B. Using our 
general scaling ansatz (13. 6p it scales for all k to [Bj, 

B nm n m ->■ (B c + a^fr~ c (d 2 + u^x)) + 0{a 2 )) ir n (3.7) 

where B c = y/r~ c (2 + s c ) is a constant operator and we have defined d = d/dx. For later 
convenience we also define 

B = d 2 + Ul (x). (3.8) 

The fact that these two terms contributing to B are of the same order fixes the scaling 
relation in (13 .6p in £ and N as functions of a. 

For a given potential the matrix H will scale to a differential operator of a fixed degree, 
say m. Assuming that H has a definite scaling dimension it will be fully determined by the 
degree m, 

H -> a< m - 2fc - 1 >/ 2 Cl c 2 (d m + ...) + o(a^ m - 2k -^ 2+1 ) = a ^-2k-i)/2 n + 0(a ( m - 2fc -i)/2 + i )) (3 Q) 

where C2 is a constant. Here we have used the result from the previous calculation that every 
power of d comes with a scaling dimension a 1//2 , and we have multiplied by the extra power 
of iV in front of the potential in (I2.30p . We may therefore write (I2.32p [B, H] = B — y as 
follows, 

G ( a (m-2fc-l)/2^ ^] + ( a (n^2fc-l)/2+l^ = fl(rf 2 + ^ _ ^ + (3 ^) 

where we have dropped the trivial constant part B c — y c that has to be satisfied. Here we 
have fixed C\ = ^/r~ c for convenience. We see that in order to get a non-trivial scaling limit 
we require m = 2k + 1. To go further let 

u = n + -{x,d} , (3.11) 
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where {, } denotes the anti-commutator. A short calculation shows that 

[B,n] = -xu[(x) +u 1 (x) - s . (3.12) 

The right hand side of (13.121) is now a multiplicative operator and we can therefore directly 
apply the technique of pseudo-differential operators reviewed in |6] to find an expression for 
% such that the left hand side of (13 . 1 2 j) is also a multiplicative operator. Indeed it is known 
that in general, 

fc+i 

H = Y,bn(B (2n - 1)/2 ) + , (3.13) 

n=l 

where here the subscript + denotes the local part of the operator fi n_1 / 2 , which contains up 
to order 2n — 1 differential operators, d 2n ~ l . The terms in (13.131) for which n < k + 1 will 
have a scaling dimension incompatible with 1-i unless the coefficients b n also come with a 
scaling dimension. The only parameter with a scaling dimension in the expression ( I2.35P is 
y and therefore the most general solution for 1-L is, 

U = c lC2 {{B^' 2 ) + - b k c 2 l s{B i2k - l ^ 2 ) + ) . (3.14) 

This form which we have deduced here on the ground of general scaling arguments is derived 
in Appendix [D] based on the relation (12.351) between H and P from the previous section. 
There it is also shown that the relative coefficient, b^ 1 , is unity. Finally we can use the 
known fact that [5] 

[ B; (^-l)/2 )+] = _ 4 l-^/ N ^ (315) 

where L\ is the Lenard differential operator and £[ denotes its d- derivative. The Lenard 
differential operator is related to the Gelfand-Dikii polynomials Ri by a trivial rescaling, 
Ci = A l R[. It is defined by the following recursion relation @, 

C' l+1 [f] = d£ l+1 [f] = (d 3 + Afd + 2f) d [/] and C = \ . (3.16) 

Because this only determines L\ one also requires the condition that £i[0] = for all /. We 
may thus write (I3.12p as 

- 4- fc Cl c 2 (£' fc+1 M - ±s£' k [ Ul \) = -xu'^x) + Ul (x) - s. (3.17) 

We now make use of the fact that under a rescaling of the variables, x — > p~ x x, u\ — > p 2 u\ 
and s — > p 2 s, we have — > p 21 Ci[u\}. Setting, 

2fc-l -1 

p = -22fe+i (cic 2 ) 2fc + 1 , (3.18) 

we obtain the equation, 

A+iK] - 4s£'Jm 1 ] = xu[{x) + 2u 1 (x) - 2s. (3.19) 

4 Note that in the notation of [6] our ui(x) — > —u(x), with B = d 2 + u\ — > Q = d 2 — u. 
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This is the differential equation our scaling function has to satisfy for the kth multicritical 
point. From [16] we see it to be a modification of the Painleve XXXIV hierarchy. For the 
first few values of k it takes the form, 

k = 1 : Ui + Quiu[ — Asu\ = xu\{x) + 2ui(x) — 2s, (3.20) 
k = 2 : uf ) + 10 (u[ (3ul + 2u'[) + mufA - 4s (uf ] + dmu'^ = xu[(x) + 2u 1 (x) - 2s, 

(3.21) 

k = 3 : u [ p + 70<u[ 3) + A2u' 1 uf ] + U Ul uf } + 70 (K) 3 + Au^u'i + u\u { ^ + 140^^ 
-4s (mS 5) + 10 (u[ (3ul + 2v![) + UtufA) = xu[(x) + 2u x {x) - 2s, (3.22) 

where we emphasise that s plays the role of a parameter and all differentials are with respect 
to x. Furthermore, from ( I2.30|) we may read off the d.s.l of P. We may write (I2.30p as, 

H = -^L(jy( B )B) + - (V'(B)B)J) - yP . (3.23) 

so y only appears in the term — yP. If we assume P scales as P — > a~ A V, where A is its 
scaling dimension we see that, 

- yP - (y c + dsa) a~ A V. (3.24) 

In order for the term containing s to match the term appearing in f!3. 141) we conclude A = 1 
and V = C2{B <y2k ^ l ^ 2 ) + . Note that the leading order term which scales as a _A must cancel 
against the other terms in (13. 23 p . If we now substitute the d.s.l of P together with the usual 
scaling ansatz into the flow equation (I2.24p and use (I3.15P we obtain, 

d s ui(x, s) = d x (x- 2C k [u x \) . (3.25) 

This equation will be useful in the next subsection in order to prove equivalence with the 
results of [12]. This equation together with (I3.19P are the main results of this subsection. 



3.2 Higher Order Analogues of Tracy Widom and Examples 

Now we consider the d.s.l of the expression (I2.8p . First we take N — > oo to obtain, 

\og[Z N (y;a,{ gi })] = N 2 I d£ (1 - log[r(£, y)\ + 0(1/N) . (3.26) 

Jo 

Note that this formula justifies setting £ c = 1; any non-analyticity of r(£,y) will only affect 
the partition function if it occurs in the range [0, 1], hence as we approach the critical point 
we must zoom in on the end point of the integration region [6j. If we now substitute the 
scaling ansatz for general k, the above takes the form, 

ro x 

logZ (fc) (s) = - / dx-{u 1 {x,s)+v 1 {x,s)) , (3.27) 

12 



where we have introduced the scaled partition function log Z^ k \s) for the fcth multicritical 
point and the limits have acquired an extra minus sign due to the rescaling (13 . 18[) . In (I3.19P 
we have only an equation for u\(x,s) and so must specify vi(x,s). Recall that we assume 
(13.51) holds for finite y and this equation depends only on the behaviour of the eigenvalue 
density at the end of the support at which the infinite wall is not imposed. It is physically 
reasonable to conjecture that since the modification due to y being finite only affects the 
local behaviour of the eigenvalue density at the other end of the support this will not have 
an effect at the opposite end and therefore (13.51) still holds. In the Appendix we show that 
the v\ does indeed satisfy (13.51) . Imposing the same boundary conditions on v\ as the y — > oo 
case, we have for the double scaling limit of the gap probability, 

f° x 

logP (fc) (s) = - / dx-(ui(x,s) -Ux(x,oo)) , (3.28) 

J -oo * 

where u\(x, oo) is the special function appearing in the expression (13 .41) i.e. it solves the kth 
member of the Painleve I hierarchy. The expression (I3.28f) together with ( 13 .191) constitute the 
main result of this paper; they give an expression for the kth multicritical gap probability and 
hence higher order analogues of the TW distribution. Indeed we can give a useful expression 
for the distribution of the largest eigenvalue using (I3.28p . we have, 

^-logPW(s) = - [ dx^d sUl (x,s). (3.29) 

dS J -oo ^ 

Using (I3T25D in (13T291 we have, 

^ logP< fc >(s) = - J° dx (c k [u x \ - ^ . (3.30) 

where we have integrated by parts assuming that the boundary conditions on u\ cause 
the boundary terms to vanish. The requirement that the boundary terms vanish implies 
that as x — > — oo we have Ui(x,s) — > Ui(x,oo) + 5ui(x,s), where 5ui(x,s) is a subleading 
contribution containing the s dependence. We note here that this boundary condition is 
compatible with (13.191) since U\(x, oo) is in fact a solution to (I3.19p . We do not discuss the 
boundary condition for 5ui(x, s) here, instead we appeal to [12] in which a proposal was 
made for the appropriate asymptotic behaviour of the solution. In the next subsection we 
relate our results to theirs; one could then in principle convert their boundary conditions 
into a condition on 5ui(x, s). 

As an example we now reproduce the standard Tracy- Widom law. For k = 1 we have 
that, 

£i[wi(x, oo)] = - m(x, oo) = | , (3.31) 
and U\(x, s) satisfies (13.191) for k = 1 which can be written as, 

«S 3) + 6^ - (x + 4sK - 2ui + 2s = , (3.32) 

where all differentials are with respect to x. Although the choice of normalisation for Ui, 
x and s in the main text gives the cleanest results for general k, it differs from the choice 
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made by Tracy and Widom [2]. We therefore perform some rescaling, in particular let 
ui(x,s) = 7 2 m(7x, 7~ 2 s) with 7 = — 2~ 1//3 . Then (13.321) becomes, 

u (3) {x, s) + 2u(x, s)(2 + 3u'(x, s)) + 2{x - 2s)u'(x, s) - As = . (3.33) 

Finally, following p3] let —2q(x, s) 2 = u(x, s) + x, substituting this into the above we find 
firstly that q(x, s) = q(x + s, 0) = q(x + s) and secondly it may be written in the form, 

q(x)W'(x) = -3q'(x)W(x) , (3.34) 

where W(x) = q"(x) — 2q(x) 3 — xq(x). This then can be integrated to show that q(x) satisfies 
Painleve II with a = 0. We now may write (I3.28P as, 

IogP«( a ) = - f dx^(u l (x,8)-Z)=-£ dx|(7 2 n(7^7^)-f)- (3.35) 
Making the change of variables x = 72 we find, 



x x 



logP^(s) = / dx — (m(x,7 2 s)+x)=— I dx xq{x + 7 2 s) 2 



logP (fc) (7 2 s) = / dx{s - x)q{x) 2 , (3.36) 

J s 

which coincides with the TW distribution up to a rescaling of s. 



3.3 Relation to Painleve II via Backlund Transformations 

In this section we relate our string equation (13 . 19[) to the Painleve II hierarchy, thereby 
making contact with the results of [T2]. We proceed using the Backlund transformation 
introduced in [16] which we now review. 

The first result we require is that if ui(x, s) satisfies the string equation ( 13.191) then, 

2C k [u 1 ] - x = 2ip{x) 2 where ^'(x)/^(x) = -W(x) , (3.37) 

and W satisfies, 

(d + 2W)C k [W' -W 2 + s] =xW + \ . (3.38) 

The proof follows from the Backlund transformation of [16], 

2C k [u x \ - x = 2^j{x) 2 and i//'(x) + {u x - s)ip(x) = . (3.39) 
Eliminating ip yields, 

K n K f ; - l -{K' n f + 2( Ul - s)K 2 n = , (3.40) 
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where we have introduced the useful quantity K n = 2C k — x. Differentiating this we obtain, 

{d 3 + Au x d + 2 Ul )C k - AsC' k = xu[ + 2u x - 2s , (3.41) 

which using the Lenard recursion formula gives (13.191) . On the other hand eliminating u% 
gives, 



2C k 



x = 2^ 2 {x) . (3.42) 



We now consider a second Backlund transformation introduced in 

2C k [W - W 2 + s] -x = 2ifj(x) 2 and ip'(x) + Wip(x) = 0. (3.43) 

Eliminating W yields again (I3.42p . whereas eliminating ip gives (13. 38 p . which is closely related 
to the Painleve II equation with a = —1/2. We can make this connection more explicit by 
using the lemma, 

c kHx ) + 4 = ' (3 ' 44) 

where z is a constant. This lemma can be verified by induction as it is shown in Appendix [El 
using the Lenard recursion relation and seems not to have appeared in the literature before. 
Hence (I3.38P can be written, 

fe-i 1 
(d + 2W)C k [W - W 2 } + t j {2^ k - 2 ^ i2k+l h) (d + 2W)C J [W - W 2 } = xW + - , (3.45) 

3=0 

where 

t-(s) = (2i + i)2( 2k - i i- 1 )/( 2k + 1 ) T{k + l/2) fc-J . 

As) ^ Zj + i)Z r(A:-j + i)r(j + 3/2) s [6Ab) 

corresponds to the Tj in [12] if k is odd, and the set of parameters U there is set to zero. In 
order to study the distribution of the largest eigenvalue the insertion of such parameters is 
not necessary, and we refer to [6] for the interpretation of these terms. 

Let us introduce the function, 

U(x,s) = -Tp(x + T ([3s),s) , (3.47) 

where j3 = 2^ k ^ 2 ^^ 2k+1 \ Using the Backlund transformations we may write the largest 
eigenvalue distribution as, 



d_ 

ds 



logPW(s) = -j"dx(c k [ Ul }- l -x^=j^ 



dx U(x — ro (/3s), s 



,2 



-MM 

dx U(x, s) 2 , (3.48) 
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and we have that U'(x, s)/U(x, s) = W(x + t (/3s), s) with W(x, s) satisfying (13.381) . Letting 
q(x, s) = —W(x + ro (/3s), s) as well as using the known symmetry of the Painleve equations 
(see e.g. [16J) we obtain the equation of [12] in q with a = +\ an d we have the relation, 

U"(x, s)/U{x, s) = q'(x, s) + q 2 (x, s) . (3.49) 

This exactly reproduces the results of [12] up to a rescaling of s. Note that we make the 
identification of our Tq(s) with the function —x(s) appearing in [12] . 



4 Conclusions and Outlook 

In this paper we have presented an alternate derivation of the higher order analogues of 
the Tracy- Widom distribution which are characterised by the spectral density vanishing as a 
(2k — l)/2-root. Instead of using Fredholm determinants we apply orthogonal polynomials to 
directly calculate a truncated partition function for the matrix model in which the upper limit 
of the integration over the eigenvalues is finite. The gap probability for the matrix model can 
then be expressed in terms of this truncated partition function. This work therefore directly 
extends an earlier paper [H] in which the usual Tracy- Widom distribution was derived. Our 
results should also be compared with those obtained earlier in [T2] via Riemann-Hilbert 
methods. We have shown the results presented here are equivalent to those appearing in [12] 
for odd values of k via a sequence of Backlund transformations. One pleasing aspect of our 
approach is that the resulting higher order analogues of Tracy- Widom can be stated more 
succinctly than in [12]. Indeed we have, 

f° x 

logP (fc) (s) = - / dx-(u 1 (x,s)-u 1 (x,oo)) (4.1) 

J—oc 2 

where p( fc )(s) is the gap probability at the kth. multicritical point for k E N, and Ui(x,s) 
satisfies, 

C' k+1 [ui] - 4:sC' k [ Ul ] = xu[ + 2u x - 2s . (4.2) 

Here all differentials are with respect to x. The issue of what boundary conditions to impose 
on U\ are discussed in the main text but are only given implicitly in terms of the boundary 
conditions chosen in [12] . 

It is worth noting that we obtain the gap probability rather than the distribution of 
the largest eigenvalue as done in [12] . One suspects then that the issue of the integration 
constant raised in [12] could be addressed more easily in our approach. 

Finally, one advantage of avoiding the use of Fredholm determinants is that one should be 
able to calculate the probability distribution for large-deviations of the maximum eigenvalue 
from its mean. Such large deviations were the main focus of [14J and since our method is a 
direct extension of theirs we expect that such an analyse could be performed for the higher 
order cases. With our approach one could also investigate issues of universality in the large 
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deviation tails. This would also give an alternative to the method in which appears to 
also be able to address these questions. This is something we hope to pursue in future work. 

Note: While completing this work an interesting paper appeared [2D] in which partial 
differential equations for Fredholm determinants associated to the d.s.l of one and two ma- 
trix models were constructed using string equations and the method of pseudo-differential 
operators. Since the gap probability may be expressed as a Fredholm determinant there is a 
large overlap in the results presented here with those appearing in (20]. However there are 
some important differences; firstly we give an explicit expression for the gap probability in 
terms of a solution to Painleve XXXIV; hence our expression is in some sense a solution to 
the PDEs derived in [20]. However, on this point, it is not immediately clear how the PDEs 
appearing in [20] compare to our expressions as they seem to require the introduction of 
extra coupling constants in the potential which also scale non-trivially in the d.s.l. Finally, 
the Backlund transformation appearing here seems to partially answer the question raised 
in j2U] concerning the relation between their approach and that of |12j . 

Acknowledgments: We acknowledge partial support by the SFB | TR12 "Symmetries 
and Universality in Mesoscopic Systems" of the German research council DFG (G.A.). 



A Equivalence of Gaussian Flow Equations 

The aim of this section is to prove the equivalence between our set of flow equations (I2.25P 
and (12.261) for the Gaussian potential and a respective set of recurrence relations that was 
derived in ref. [H] in eqs. (42) and (41) there: 

Sn-Sn-i = -i<Vog[r n ] , (A.l) 
la 

r n+ i - r n _i + s 2 n - s 2 n _ x = -d & log[r n ) , (A. 2) 



with a = Na. Obviously the first equation agrees with (I2.25p . To show that eq. flA.2j) also 



follows from our formalism consider first the matrix P in the form of eq. (I2.19p . Due to its 
antisymmetry the diagonal has to vanish, 

P nn = A nn + C nn - aB nn = --d y log[h n ] - as n = 

=> s n = e- & y\ n {yf , (A.3) 

2a 

where we have used (12.121) and (I2.17P or equivalently (I2.2ip in the second line. Our main 
relation is obtained from the lower diagonal, 

P n , n -i = Ai,n-i + Cn,n-i - «S n ,n-i = - e~ ay iT n (y)n n -i(y) - ay/r^ 

-Pri.— l.ri. \l fn. 



n = j7~ n e ay 7r n (y)7r n _i(y) + 2ar n , (A.4) 
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using eqs. (I2.13P and (I2.2ip as well as the antisymmetry of P. We can now use the three 
step recurrence relation to eliminate y/r^7T n -x(y) as well as use ( IA.3j) : 



n 



<^> 2n + 1 



-ay 2 



e n n (y) {{y - s n )n n {y) - y/r n+1 ir n+1 (y)) + 2ar n 
-2a(y - s n )s n - (n + 1 - 2dr n+1 ) + 2dr n , 
2a((s n - y)s n + r n+1 + r n ) , 



(A.5) 



where we have also applied flA.41) for n — > n + 1 . Taking the difference between eq. flA.5j) at 
rz and n-1 we finally arrive at 



2d(y(s n _i - s n ) + s 2 n - 4_i + r n+1 - r. 



n-li 



2a{y—d y log[r n ] + 4 ~ s n-i + r «+i - r n- 
2a 



(A.6) 



upon inserting eq. (lA.lj) . The last step to arrive ateq. fTA~2|) requires the following identity: 



-ad & log[r r , 

It follows by comparing the y- and d- derivative of 



1 



1 



dz e 



1 



dX e~ 7r n (A) 2 

) 

For the former we have already from (1A.3|) 

^log[/i„]= e-^ 2 7r„(y) 2 , 
whereas for the latter we obtain 

= ~ + #r e-^vr^y) 2 - ^log[^ 
2a 2a 



n-l- 



i 

'2d 



(A.7) 



(A. 



(A.9) 



(A.10) 



Inserting these equations into each other and taking the difference log[/i n ]— log[/i n _i] = log[r r , 
we arrive at (IA.7[) . 



B Multiplying Upper and Lower Triangular Matrices 

Any matrix D can be decomposed into its strictly upper (+), diagonal (d), and strictly lower 
triangular (— ) part: 

D = D + + D d + D_ (B.l) 

The aim of this appendix is to express the + and — part of the product of two matrices in 
terms of the +, d and — parts of the individual factors, in order to express matrix H in 
terms of matrix P in the main text eq. (12. 35 p . Obviously multiplying any matrix by the 
diagonal part of a matrix D d will not change the character of that matrix to be +, d, or — : 

(BD d ) + = B + D d , (BD d ) d = B d D d , (BD d )_ = B_D d . (B.2) 
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Furthermore multiplying two strictly upper (lower) triangular matrices gives a strictly upper 
(lower) triangular matrix, (B + D + ) + = B + D + and (B-D-)- = B-D-. 

From now on we will use that our matrix B eq. (12.1 Op is tridiagonal. Therefore it can 
only shift off-diagonals at most up or down by one, that is 

{BD + ) + = B + D + + B d D + + (B_D + ) + = BD + -(B_D + ) d , (B.3) 

due to (BD + ) d = (B^D + ) d and (£>_.D + )_ = 0. Analogously one can deduce 

(££>_)_ = (B+£L)_ + B d D_ + = BD_ - {B + D_) d , (B.4) 

because of (BD_) d = (B + DJ) d and (B + DJ) + = 0. Consequently we can write down for 
any product BD its respective +, d and — part: 

(BD) + = BD + -(B_D + ) d + B + D d , 
(BD) d = (B + D_) d + B d D d + {B_D + ) dl 

(££>)_ = BD_ - (B + D_) d + B_D d . (B.5) 

We can now proceed expressing H through P. Defining 

Na 

F = --V'(B) (B.6) 

which commutes with B and is symmetric as is B, we have 

P = F + -F_, (B.7) 

which has no diagonal part, P d = 0. Finally we obtain 

H = (BF) + -(BF)_-y(F + -F_) 

= (B - y)(F+ - F_) - (5_F + ) d + (£?+F_) d + B + F d - B_F d 

= {B-y)P-{BP) d + {B + -B_)F dl (B.8) 

where in the last step we have added trivially vanishing terms (B_PJ) d = = (B + P + ) d and 
(B d P) d = in order to express everything in terms of P. 



C The Double Scaling Limit in the Case y — > oo 

Here we recall some standard results about the d.s.l of the Hermitian one-matrix model 
whose partition function is Z N (oo; a, {g{\). The equations (12.61) and (12. 9p hold in this case 
with now all quantities independent of y. The recursion coefficients in the definition of (12. 9p 
are determined by the recursion relations, also known as the string equations, 

aV'(B) n , n+1 = n/N , 

V'{B) nn = 0. (C.l) 
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Note that these equations can be obtained from (I2.19P by setting derivatives of y, and hence 
C, to zero and using (12.131) . By taking the large N limit we can rewrite these equations in 
terms of the continuous variable £ = n/N. It is known (see e.g. [7]) that for certain critical 
potentials in which the eigenvalue density vanishes as a (2k + l)/2-root at the right and a 
(2/c' + l)/2-root at the left end of the support, respectively, with k' < k, that substituting the 
ansatz (13. ip and (13.21) into fIC.ll) results, after the rescaling (13.181) . in the equations (13.41) and 
(13.51) . It was then shown in [8] that the correct solution to (13. 5p was v\ = 0. The equations 
(IC.ip may also be written as, 



[B,P] = 1, (C.2) 

where P is defined by (I2.22p . The equations (13 .4p and (I3.5P arise in this operator formalism 
since it is easy to show that B scales to a differential operator independent of k, 

B B c + aci (d 2 + ui{x)) + 0(a 2 ) =B C + ac x B + 0(a 2 ) . (C.3) 

Furthermore, due to the finite number of nonzero off-diagonals in P and its antisymmetry, 
we expect P to scale to an anti-hermitian differential operator. In general for k' < k, 
P — > a~ v V + . . . where rj > 1. Requiring that we have a scaling rj = 1 for the leading 
term determines the functions Wj, in particular it gives (I3.5p . Given that [B, V] must be a 
multiplicative operator we have V oc (B 2k ~ l ) + and therefore (13 .4p . up to rescalings. 



D Derivation of the d.s.l for H 

In contrast to the main text in subsection 13.11 we first derive here the scaling of operator 
P, and then using the relation between H and P show the d.s.l for H to be of the form 
(I3.14p . We now consider the case of the potential in Z N (y; a, {gi}) taking the form of a 
critical potential of Z^(oo; a, {gi})- The left hand side of the flow equation (I2.24p is exactly 
the left hand side of (1C 2|) and the right hand side of (12.241) has the same scaling dimension 
as the right hand side of (1C.2I) . This means that using the ansatz (I3.6P and (13.21) . we have 
for k' < k that v i satisfies (13.51) and furthermore we can also conclude 

P -> a- l V + 0(1) = a- 1 c 2 (B {2k - 1)/2 ) + + 0(1) . (D.l) 

The double scaling limit of the flow equation is then, using (13.151) 

l-d a u 1 = -A 1 - k c 1 c 3 £! k [u 1 ] , (D.2) 

which will be of use in the following. We now turn our attention to H in (12.281) and prove, 

U = c 2Cl {(B^' 2 ) + - s(B^/ 2 ) + ) + i{x, d} . (D.3) 
Recall the expression (I2.35p . we have, 

H=±{B,P}-yP-±{B + -B^\V'{B)] A } . (D.4) 
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If we now consider the d.s.l it is a simple matter to prove that independent of k, 

B + - B_ -> -2 Cl a 1/2 d + 0(a). (D.5) 

To find the double scaling limit of the operator ^[V (B)]^ we must work slightly harder. 
Recalling the expression for P f 1 2 . X 9 j) and (I2.17p . then the antisymmetry of P implies, 

/ — Na 

C nn = -d v log[y/K] = —V'(B) nn , (D.6) 

which gives, 

- l -d y log[r n ] = ^V'(B) nn - ^V'(fl)n-i,»-i . (D.7) 

where V'(B) nn are trivially related to the operator [V'(.B)]d- In the large N limit we have 
?Y L V'(B) nn — > V(£,y), where V(£,y) is a continuous function of £. We therefore have, 

-\a K ]og[r(t,v)) = jjd ( fr(t,y). (D.8) 

Substituting in (13 .6p and (13 .2p while assuming an arbitrary scaling dimension for V i.e. 

V(£, y) a"V(a;, s) gives, 

^-<9 S s) + ^(s, s)) = ^d sUl & a) = a^ +1 / 2 d x V(x, s), (D.9) 

from which we see see that rj = —1/2. Furthermore, using (ID.2[) we have, 

V(x, s ) = J-(x + ^^dCaAK]), (D.10) 
4ci 

where we have integrated and assumed the integration constant is zero. 

Returning now to f]D.4[) . the form of the d.s.l for all operators is now known. Substituting 
them in we obtain, 

H = ^c lC2 {B, {B k - l ' 2 ) + } - c lC2 s{B {2k - l ^ 2 ) + + ^{x, d} + A- k Cl c 2 {C k , d}, (D.ll) 
which upon using the identity [Sj {B k+1 l 2 ) + = \{B, (£ fe - 1/2 ) + } + 4r k {C k , d} becomes fHX3l) . 



E Shift Identity for Lenard Differential Operators 

We claim that when z is a constant it holds for any integer k G N: 

C k \u(x) + z} = V (4z) k -i— F(fc+ X 1 { 2) —CAuix)] . 

L f^ Q T(k-j + l)F(j + 1/2) n v n 
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We now prove this via induction. Firstly, it is trivially true for k = 1. Now let us assume, 

k 

C k [u(x) + z] = Y,{^) k ~ j o i f ) C j [u{x)) , (E.2) 

where we have defined 

J T{k - j + l)r(j + 1/2) 

and for any fcGE It is simple to see that the following identity holds for these constants: 

af +1) =af + afl 1: for j = l,...,k . (E.4) 



Let us now substitute the assumption (1E.2I) at k + 1 into Lenard's recursion relation (13.161) : 

fc+1 
fc+i 

+ = ^2af +1 \iz) k+1 - j C j [u)+C k+1 . (E.5) 
j'=i 

Here we have introduced C k+ \ as an integration constant and used the fact that Cq is a 
constant as well as that a k k ^ = ol£' = 1. To complete the proof by induction we must fix 
Ck + i. This can be done by using the fact that Cj[0] = for all j > 0, so that we have, 

C k+1 [z) = C k+1 . (E.6) 

Due to the homogenous scaling of we know there exists only one term in 

which is non-zero when u(x) = z = const; it is, f3 k+ iu k+1 , where (3 k+ i is a known constant 

that follows from the properties of the Lenard differential. We therefore have, 

C k+l = (3 k+1 z k+l . (E.7) 

Let us finally determine k . This can be done by noting that in the Lenard recursion relation, 
the only term in C k which contributes to the u k+l term in £ k+ i, is (3 k u k . Hence we have 
from the Lenard recursion relation, 

d((3 k+1 u k+1 ) = AuP k u k - l u' + 2u'p k u k = A^^p k d{u k+l ), (E.8) 

* Pk+1 ~ 4 r(^ + 2)r(3/2) /31 ' (K9) 

where j5\ = 1. Substituting the known value of C k+ \ into our previous relations using 
C = 1/2 gives, 

fc+i 

C k+1 [u + z] = ^«f +1) (4^) fc+1 -^-[n] , (E.10) 

j=0 

thereby completing the proof. 
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